Applying Streaming Analytics to Tomographic Imaging

• Imaging in 2015
• Traditional Approaches
• Streaming Analytics
• Drilling Down
• Outlook / Developments

X-Ray

• Swiss Light Source (SRXTM) images at (>1000fps) \( \rightarrow \) 8GB/s
• Diffraction patterns (cSAXS) at 30GB/s
• Nanoscopium (Soleil), 10TB/day, 10-500GB file sizes, very heterogenous data
• Swiss Radiologists will make 1PB of data year

Geospatial

• New satellite projects (Skybox, etc) will measure >10 petabytes of images a year

Personal

• fps1000 - (640 x 480 x 840 fps) for /$400 \( \rightarrow \) 780Mb/s \( \rightarrow \) 66TB/day

  • GoPro 4 Black - 3840 x 2160 x 30fps for $600 \( \rightarrow \) 750Mb/s \( \rightarrow \) 64TB/day \( \rightarrow \) 23PB/year

The breakdown of science into component tasks has changed from a time perspective over the last years and how we expect it to change by 2020.

• Experiment Design (Pink)
• Data Management (Blue)
• Measurement Time (Green)
• Post Processing (Purple)

The primary trends to focus on are the rapid decrease in measurement time and increasing importantance of post-processing.

The breakdown of science into component tasks has changed from a time perspective over the last years and how we expect it to change by 2020.

• Experiment Design (Pink)
• Data Management (Blue)
• Measurement Time (Green)
• Post Processing (Purple)

The primary trends to focus on are the rapid decrease in measurement time and increasing importantance of post-processing.

• Measure a full mouse brain (1cm\( ^3 \)) with cellular resolution (1 \( \mu \) m)
• 10 x 10 x 10 scans at 2560 x 2560 x 2160 \( \rightarrow \) 14 TVoxels
• 0.000004% of the entire dataset

• 14TVoxels = 56TB
• Each scan needs to be registered and aligned together
• There are no computers with 56TB of memory
• Analysis of the stitched data is also of interest (segmentation, vessel analysis, distribution and network connectivity)

Thank you Alessandra and Alberto!

The first of challenges with these data is to start stitching them together.

• Algorithmically it is very simple:

• \( \textrm{Overlap}_{A \rightarrow B} = \)

• \( \textrm{argmax}_{\vec{x}}\left[\mathcal{FT}^{-1} \left\{\mathcal{FT}\{A\}\mathcal{FT}\{B\}^{*}\right\} \right] \)

Computationally

We need to stitch each scan with all of its neighbors

• \( x=\pm 1 \)
• \( y=\pm 1 \)
• \( z=\pm 1 \)

We have a few more parameters (\( \theta \)) and we know the initial positions so we optimize in a limited range rather than the entire image.

• \( \rightarrow \) we need to add another loop
• Iterate as long as the offset vectors are changing

Now that our task is clear we just execute \( \rightarrow \)

Super workstation

• 160GB RAM
• awesome-sauce graphics card
• 40 cores with hyperthreading
• SSD drives

Matlab + FIJI

• Premade algorithms
• Just click and go

• 10x10x10 scans x 8 neighbors \( \rightarrow \) 8,000 registrations per iteration
• 20 - 40 minutes per registration
• 1000-4000 hours \( \rightarrow \) 40-160 days per iteration

Basic Topology

• Driver Node

The machine being used by the user which is responsible for creating jobs which need to run

• Master Node

Distributes task across multiple machines, coordinates communication between machines

• Worker Nodes

The nodes where the compuation is actually done

• \( \rightarrow \) Tell the driver to load all images
• The driver determines which data need to be loaded
• Driver \( \rightarrow \) Master exact data to load

• tell machine 1 to load images 0 to 2
• tell machine 2 to load images 3 to 5
• …

tell share the images to calculate the overlap

• \( \rightarrow \) tell machine 1 to send image 0 to machine 2
• \( \rightarrow \) tell machine 2 to calculate: \( \textrm{Overlap}(A\rightarrow B) \)
• …

• Every task on every machine must be explicitly specified
• Substantial effort in coordination
• Huge amounts of data traffic between shared storage and machines
• Not flexible or elastic

Correctness

The most important job for any piece of analysis is to be correct.

• A powerful testing framework is essential
• Avoid repetition of code which leads to inconsistencies
• Use compilers to find mistakes rather than users

Easily understood, changed, and used

Almost all image processing tasks require a number of people to evaluate and implement them and are almost always moving targets

• Flexible, modular structure that enables replacing specific pieces

Fast

The last of the major priorities is speed which covers both scalability, raw performance, and development time.

• Long waits for processing discourages exploration
• Manual access to data on separeate disks is a huge speed barrier
• Real-time image processing requires millisecond latencies
• Implementing new ideas can be done quickly

• Untested Algorithms

  • 5 different implementations (GPU, cluster-version, local version, Matlab prototype version)

• Image Definitions

  • Too Strict
  • Mat: CV_8UC1 1D array of unsigned bytes
  • Img<RealType< T > & NativeType< T > >> 2D array of generically typed with cursor access
  • itk::Image< TPixel, VImageDimension >
  • Too Loose
  • imread could be 2D, maybe 3D, sometimes even 4D
  • numpy runtime problem determination

• Many independent components
• Complicated workflows
• Fragile glue

• Courtesy of
  • https://d13yacurqjgara.cloudfront.net/users/13496/screenshots/1333062/attachments/188364/analytics-app-2-animated.gif
  • https://g.twimg.com/blog/blog/image/TAD_Buzzfeed_Animated.gif

• You can analyze streaming tweets during World Cup matches to see goals before you see them on television

• You can process a petabyte of clicks to optimize ad-campaigns

• Powerful tools

• Distributed computation

• Interactive computation

Boring problems

Velocity, Volume, Variety

When a ton of heterogeneous is coming in fast.

Performant, scalable, and flexible

When scaling isn't scary

10X, 100X, 1000X is the same amount of effort

When you are starving for enough data

Director of AMPLab said their rate limiting factor is always enough interesting data

O 'clicks' per sample

Question not implementation driven

• Ask the questions you should not the questions you can

You conduct the orchestra, you don't spend time telling each person what to do

Big Data is a very popular term today. There is a ton of hype and seemingly money is being thrown at any project which includes these words. First to get over the biggest misconception big data isn't about how much data you have, it is how you work with it

Before we can think about Big Data we need to define what it's replacing, Small Data.

Small Data

The 0815 approach, using standard tools and lots of clicking, one image at a time

Storage

• 8TB for 2 years is \( \approx \) $5000 (> petabytes)
• From SLS to Cloud \( \rightarrow \) 231MB/s
• Standard HD would cost $1200
• K. Mader, et al, Moving image analysis to the cloud: A case study with genome-scale tomographic data, AIP Conference Proceedings, X-Ray Microscopy

Processing Power

Elastic, On-Demand

• Rather than building an analysis as quickly as possible and then trying to hack it to scale up to large datasets
  • chose the framework first
  • then start making the necessary tools.
• Google, Amazon, Yahoo, and many other companies have made huge in-roads into these problems
• The real need is a fast, flexible framework for robustly, scalably performing complicated analyses, a sort of Excel for big imaging data.

Apache Spark

• World Record Holder for Terabyte and Petabyte Sort
  • 100TB of data on 206 machines in 23 minutes
  • 1PB of data in 190 machines in <4 hours
• Scalable
  • >10000 computers
  • storing and processing exabytes of data
  • sustained data rates of 3Gb/s/node (map) and 1.1Gb/s/node (reduce)
• fault tolerance
  • 2/3rds of computers can crash and a request still accurately finishes
• hardware and software platform indpendence (Mac, Windows, Linux)

These frameworks are really cool and Spark has a big vocabulary, but flatMap, filter, aggregate, join, groupBy, and fold still do not sound like anything I want to do to an image.

I want to

• filter out noise, segment, choose regions of interest
• contour, component label
• measure, count, and analyze
• …

Spark Image Layer

• Developed at 4Quant, ETH Zurich, and Paul Scherrer Institut
• The Spark Image Layer is a Domain Specific Language for Microscopy for Spark.
• It converts common imaging tasks into coarse-grained Spark operations

We have developed a number of commands for SIL handling standard image processing tasks

Fully exensible with

New detectors certainly provide a substantial amount of new data, but what about old data?

• At PSI we have years of tape archives with thousands of interesting experiments.
• At hospitals and satellite imaging companies it is even more

PI:

Look at this new technique from famous journal, let's try it on all of our old data

PhD Student:

1. Find old data
2. Copy old data one disk at a time to current machine
3. Read in old data, Run analysis, Save results
4. Repeat from 2.

2 years later

It seems to give the same results

PI

Well what about this technique?

Repeat ad naseaum

How many of the bone samples measured had a cortical thickness larger than 30mm? \[ \downarrow \]

SELECT COUNT(*) FROM (
  SELECT THICKNESS(
    SEGMENT_CORTICAL(
      GAUSSFILT(imgData)
    )
  ) AS thk WHERE thk>30
)

Is there a correlation between cortical thickness and mineral content?

\[ \downarrow \]

SELECT THICKNESS(cortical),
  IMG_AVG(cortical) FROM (
  SELECT SEGMENT_CORTICAL(
    GAUSSFILT(imgData)
    ) AS cortical
  )

The query is processed by the SQL parsing engine

• sent to the database where it can be combined with locally available private data
• sent to the image layer for all image processing components to be integrated
• sent to the Spark Core so specific map, reduce, fold operations can be generated
• sent to the master node where it can be distributed across the entire cluster / cloud

In the biological imaging community, the open source tools of ImageJ2 and Fiji are widely accepted and have a large number of readily available plugins and tools.

We can integrate the functionality directly into Spark and perform operations on much larger datasets than a single machine could have in memory. Additionally these analyses can be performed on streaming data.

val wr = new WebcamReceiver()
val ssc = sc.toStreaming(strTime)
val imgList = ssc.receiverStream(wr)

Filter images

val filtImgs = allImgs.
  mapValues(_.run("Median...","radius=3"))

val totImgs = inImages.count()
val bgImage = inImages.
  reduce(_ add _).multiply(1.0/totImgs)

val eventImages = filtImgs.
    transform{
    inImages =>
      ...
      val corImage = inImages.map {
        case (inTime,inImage) =>
          val corImage = inImage.subtract(bgImage)

val corImage = inImages.map {
      case (inTime,inImage) =>
        val corImage = inImage.subtract(bgImage)
        (corImage.getImageStatistics().mean,
          (inTime,corImage))
    }

val eventImages = filtImgs.
    transform{
    inImages =>
      val corImage = inImages.map {
        case (inTime,inImage) =>
          val corImage = inImage.subtract(bgImage)
          (corImage.getImageStatistics().mean,
            (inTime,corImage))
      }
      corImage
  }

eventImages.filter(iv => Math.abs(iv._1)>20).
  foreachRDD(showResultsStr("outlier",_))

A spinoff - 4Quant: Big Image Analytics

• Bringing Big Data and Science Together
• Quantitative Search Machine for Images
  • Make analyzing huge sets of image data like a database search
• Custom Analysis Solutions
  • Custom-tailored software to solve your problems
  • Streaming / Acquisition Solutions
  • Integration with Image Database Systems
  • Cloud Solutions

Education / Training

• Consulting
  • Advice on imaging techniques, analysis possibilities
• Education
  • Workshops on Image Analysis
  • Courses / Training
  • Quantitative Big Imaging
• Contact Us: info@4quant.com

While our framework is commercial, we build on top and integrate into open-source tools so that the research is reproducible by anyone, anywhere, using any machine (it just might take a bit longer)

• ImageJ2
• KNIME
• BigDataViewer

We also try to publish as much code, educational material, and samples as possible on our github account.

• https://github.com/4Quant

Big Data Finite Element Analysis using GraphX and Spark

• https://github.com/4Quant/spark-fem

• https://github.com/kmader/Quantitative-Big-Imaging-2015/blob/master/Exercises/10-Description.md

• Built on Freely-available, Open Source tools
• Validation with synthetic datasets
• Understanding and improving performance

Here we use a KNIME -based workflow and our Spark Imaging Layer extensions to create a workflow without any Scala or programming knowledge and with an easily visible flow from one block to the next without any performance overhead of using other tools.

• AIT at PSI and Scientific Computer at ETH
• TOMCAT Group

We are interested in partnerships and collaborations

Learn more at

• 4Quant: From Images to Statistics - http://www.4quant.com
• Quantitative Big Imaging Course at ETH Zurich
• X-Ray Imaging Group at ETH Zurich - http://bit.ly/1gD8wKb

• First we start our cluster:
  spark-ec2.py -s 50 launch 4quant-image-cluster

• Load all of the samples (56TB of data)

loadImage(
  "s3://../brain_*_*_*/rec.tif"
  )

\[ \downarrow \] A Resilient Distributed Dataset (RDD) \[ \textrm{Images}: \textrm{RDD}[((x,y,z),Img[Double])] =\\ \left[(\vec{x},\textrm{Img}),\cdots\right] \]

• Now start processing, run a gaussian filter on all the images to reduce the noise

filteredImages = Images.gaussianFilter(1,0.5)

• Calculate the volume fraction at a threshold of 50%

volFraction = Images.threshold(0.5).
  map{keyImg =>
    (sum(keyImg.img),keyImg.size) 
    }.reduce(_+_)

Developing in Scala brings additional flexibility through types[1], with microscopy the standard formats are 2-, 3- and even 4- or more dimensional arrays or matrices which can be iterated through quickly using CPU and GPU code. While still possible in Scala, there is a great deal more flexibility for data types allowing anything to be stored as an image and then processed as long as basic functions make sense.

[1] Fighting Bit Rot with Types (Experience Report: Scala Collections), M Odersky, FSTTCS 2009, December 2009

What is an image?

A collection of positions and values, maybe more (not an array of double). Arrays are efficient for storing in computer memory, but often a poor way of expressing scientific ideas and analyses.

• Filter Noise?

combine information from nearby pixels

• Find objects

determine groups of pixels which are very similar to desired result

We have all of the filtered images and we want to stitch them together in a smart way.

• Start simple and match the images up with its neighbors (all image permutations and filter out the proximal ones)

pairImages = Images.
  cartesian(Images).
  filter((im1,im2) => dist(im1.pos,im2.pos)<1)

The cross correlation can then be executed on each pair of images from the new RDD (pairImages) by using the map command

displacementField = pairImages.
  map{
  ((posA,ImgA), (posB,ImgB)) =>

    xcorr(ImgA,ImgB,
      in=posB-posA)

  }

If you looked at one 1000 x 1000 sized image every second

It would take you
139 hours to browse through a terabyte of data.

The standard approach to producing data is top-down.

• Selecting a model
• generating hypotheses
• conducting experiments
• updating the model.

In a study with 80 samples with a variability of 10% and a 5% difference between them, we can afford only 1 test and maintain the commonly accepted \( p<0.05 \)

We thus need smart people making and choosing the best models so that we can preserve the signifance.

• Rather than impose a preconceived, possibly biased, model on data we should investigate what possible models, interpretations, or correlations are in the data (possibly in the phenomena) that might help us understand it. - Michael Brodie: MIT and Former Chief Scientist at Verizon

The Big Data era suggests a new approach. First given an abundance of data, we can afford to run many tests.

In a study with 800 samples with a variability of 10% and a 5% difference between them, we can afford over 10M tests and maintain the commonly accepted \( p<0.05 \)

For knowledge creation this means a switch

• from long, arduous studies involving careful model selection
• to hypothesis and model creation based on the correlations and trends found in the raw data

Moores Law

\[ \textrm{Transistors} \propto 2^{T/(\textrm{18 months})} \]

Based on data from https://gist.github.com/humberto-ortiz/de4b3a621602b78bf90d

There are now many more transistors inside a single computer but the processing speed hasn't increased. How can this be?

• Multiple Core
  • Many machines have multiple cores for each processor which can perform tasks independently
• Multiple CPUs
  • More than one chip is commonly present
• New modalities
  • GPUs provide many cores which operate at slow speed

Parallel Code is important

The range of cloud costs (determined by peak usage) compared to a local workstation with utilization shown as the average number of hours the computer is used each week.

• People usually only work ~40 hours a week (in europe at least)
• Adding IT costs and management easily shifts the blue cruve substantially upward.

• ImageJ macro for segmentation
• Python script for shape analysis
• Paraview macro for network and connectivity
• Python/MPI script to pool results
• MySQL Database storing results

1.5 man-years

Advantages

• Faster than manual methods
• Less bias
• More reproducible
• Can be highly optimized / performant

Disadvantages

• Compartmentalized analysis
• Time-consuming
• Complex interdependent workflow
• Fragile (machine crashes, job failures)
• Expensive - Measure twice, cut once, not exploratory

For many traditional scientific studies, the sample count is fairly low 10-100, for such low counts, if we look at the costs (estimated from our experience with cortical bone microstructure and rheology)

• Development costs are very high for Big
• Traditionally students are much cheaper than the 150kCHF / year used in this model
• Medium Data is a good compromise

As studies get bigger (>100 samples), we see that these costs begin to radically shift.

• The high development costs in Medium and Big \( \rightarrow \) 0
• The high human costs of Small and Medium \( \propto \) \( Samples \)

Bottleneck is filesystem connection, many nodes (10+) reading in parallel brings even GPFS-based infiniband system to a crawl

One of the central tenants of MapReduce™ is data-centric computation \( \rightarrow \) instead of data to computation, move the computation to the data.

• Use fast local storage for storing everything redundantly \( \rightarrow \) less transfer and fault-tolerance
• Largest file size: 512 yottabytes, Yahoo has 14 petabyte filesystem in use

Combining many different components together inside of the Spark Shell, IPython or Zeppelin, make it easier to assemble workflows

Web and Image Database Examples